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Abstract 

The approximate representation of a quantum sohd as an equivalent composite semi-classical 
solid is considered for insulating materials. The composite is comprised of point ions moving on 
a potential energy surface. In the classical bulk domain this potential energy is represented by 
pair potentials constructed to give the same structure and elastic properties as the underlying 
quantum solid. In a small local quantum domain the potential is determined from a detailed 
quantum calculation of the electronic structure. The formulation of this description as a sequence 
of physical approximations is considered in some detail. The primary new ingredients are 1) 
a determination of the pair potential from quantum chemical data for equilibrium and strained 
structures, 2) development of 'pseudo-atoms' for a realistic treatment of charge densities where 
bonds have been broken to define the quantum domain, and 3) inclusion of polarization effects on 
the quantum domain due to its environment. This formal structure is illustrated in detail for an 
Si02 nanorod. For each configuration considered, the charge density of the entire solid is calculated 
quantum mechanically to provide the reference by which to judge the accuracy of the modeling. 
The construction of the pair potential, the rod, the pseudo-atoms, and the multipoles is discussed 
and tested in detail. It is then shown that the quantum rod, the rod constructed from the classical 
pair potentials, and the composite classical/quantum rod all have the same equilibrium structure 
and response to elastic strain. In more detail, the charge density and forces in the quantum 
subdomain are accurately reproduced by the proposed modeling of the environmental effects even 
for strains beyond the linear domain. The accuracy of the modeling is shown to apply for two quite 
different quantum chemical methods for the underlying quantum mechanics: transfer Hamiltonian 
and density functional methods. 
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I. INTRODUCTION 



The field of multiscale modeling has opened a new era to computational science for study- 
ing complex phenomena like fracture, hydrolysis, enzymatic reactions, solute-solvent studies, 
hydrogen embrittlement, and many other chemo-mechanical processes in macroscopic sam- 
ples. Often these phenomena require a very accurate description at one scale, while at 
another scale a much coarser method can be applied to get satisfactory results. In fact, the 
coarser description for the bulk sample is required because the more fundamental methods 
are computationally too intensive to be applied to the entire system. These scales may be 
length or time, or a combination of both. This scheme of combining different models at 
different scales to achieve a balance of accuracy, efficiency, and realistic description is known 
as multiscale modeling. It is accomplished by applying the high accuracy method only in 
a small domain where it is needed the most and more approximate methods for the rest of 
the bulk where they are appropriate. For example, in crack propagation one has to apply 
detailed quantum mechanics at the tip of the crack. Here the bonds are breaking between 
the atoms, causing a marked deformation of the electron cloud and charge transfer between 
ions. But far from the crack tip where the atoms are less deformed they can be described by 
classical mechanics with appropriate potentials. Thus, the problem arises of linking quite 
different descriptions for the different length scales^. 

There is a broad diversity of other examples requiring multiscale modeling: chemo- 
mechanical polishing^, tidal wave prediction^, atmospheric sciences, embrittlement of nu- 
clear reactors^, and many biological systems^. Two classes can be distinguished^, "serial 
multiscale modeling" where the various scales are weakly coupled but the computation of 
parameters at smaller scales is required for its use in more phenomenological models at a 
larger scale, and "concurrent multiscale modeling" where the various descriptions applied 
on different scales should all be nested with proper boundary conditions. The multiscale 
modeling discussed here is of the "concurrent" type, although there are components that 
can be considered as the "serial" type (e.g., parameterization of the underlying approximate 
quantum method used). 

Three different length scale levels are distinguished, the nanoscale where the details of 
electron structure and quantum chemistry are important, the atomic or microscale where ap- 
propriate classical pair potentials capture the structure accurately, and the macroscale where 
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a continuum mechanical description applies. The present work addresses only one partic- 
ular subclass of multiscale modeling, known as hybrid quantum mechanical (QM) /classical 
mechanical (CM) problems, and their application to solid insulators - silica in particular. 
The further embedding of the atomic scale in a macroscale model is not considered here. 
Early work on this QM / CM topic began with Warshel and Levitt^- and has accelerated over 
the past decade^. Applications of QM/CM include biological systems where this scheme 
is sometimes referred to as QM/MM (molecular mechanics) modeling (enzymes, DNAs 
and proteins)^'i2iiiii24H^ vibrational spectroscopy^^, electronic excitations^ i^^i^^ , hydrolysis 
of silic a^^i^^ , and solute-solvent problemsi^i^Siii. 

The challenge of the QM/CM simulation is to move from one length scale to another, 
as smoothly as possible. A solid or large cluster behaves as a single "molecule" so the 
partitioning becomes more complicated due to covalent bonds that are cut between the 
QM and CM regions. These dangling bonds give rise to incorrect charge densities in the 
QM domain and other pathologiest^SiS^. A multiscale modeling scheme is proposed here 
that addresses this and other problems of the quantum embedding to provide a means for 
faithful coupling between the QM and CM regions. There are many tests for fidelity across 
the QM/CM interface, such as bond angles, bond lengths, and proton affinities. Instead, 
the criteria here are preservation of accurate electronic charge densities and total forces in 
the QM sub domain. These are the physically significant properties based more closely in 
the theoretical structure being modeled, as described in the next section. The approach 
includes the following two components: 

(i) Modeling environmental effects for the quantum domain in an accurate and simple 
manner. Here, accurate means that the forces and charge densities of the QM domain 
are reproduced to within a few percent. This entails a chemically correct treatment of 
the dangling bonds, and an account of longer range Coulomb influences from the classical 
domain. Simplicity, means that the approximations made are transparent and that the 
difficulty of the quantum calculation is not increased. Here, the valency problem is solved 
by pseudo-atoms and the longer range environmental effects are included through low order 
polarization effects^i. 

(ii) Developing a new classical potential for the classical region tailored to the properties 
of interest. A new potential for the entire CM region is designed, based on selected force data 
calculated from the quantum method used in the QM domain. This data is generated for 
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both equilibrium and near equilibrium states. The test for a suitable potential, in addition 
to producing the right structure, is the accurate reproduction of the chosen property to be 
studied, for near equilibrium states (typically some linear response characteristic)^. The use 
of this new potential prevents any mismatch of that response property across the QM/CM 
boundary. Hence different problems (fracture or corrosion) might require the construction 
of different potentials for the same system. 

The criterion for the success of (i) and (ii) is that calculations of the desired properties 
should be indistinguishable for the quantum system, purely classical system, and the com- 
posite QM/CM system for near equilibrium states. Only then are the desired applications 
to far from equilibrium states reliable. 

In the next four sections, the theoretical basis for this modeling scheme is described in 
some detail. The object is to provide a means to assess the approximations made more 
directly, and also to provide the basis for their extension to other systems. The complete 
formalism is then tested for a model system - a silica nanorod^^ (Figure - in section VI. 
This silica nanorod containing 108 nuclei has been chosen since a quantum treatment of the 
whole sample is possible to assess quantitatively all aspects of the proposed scheme. Still, 
the system is big enough to yield bulk properties like stress-strain response (also shown in 
Figure P). The nanorod has the proper stoichiometric ratio of silicon to oxygen observed in 
real silica (1:2) and is considered a viable model for studying silica. One of the rings near 
the center of the rod is chosen to be the QM domain and rest of the rod is treated classically 
(see Figure El below) . 

The system is studied under uniaxial strain. The quantum mechanical method used is the 
Transfer Hamiltonian (TH) method^ using an NDDO type theory parameterized by coupled 
cluster data (described more completely in section VI). The construction of the classical pair 
potential is described and shown to yield the same Young's constant as obtained from the 
QM calculation, to within a few percent. Next, the training of the pseudo atoms and 
validity of polarization representation for the environment of the QM system is tested. The 
resulting force and charge densities for the QM domain are accurately reproduced to within 
a few percent. Finally, the composite QM/CM rod is shown to have the same structure and 
elastic properties as both the QM and CM rods for near equilibrium states, as required. The 
entire analysis is repeated using density functional theory (DFT) as the underlying quantum 
theory, instead of the TH quantum mechanics. The same level of accuracy for the modeling 
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FIG. 1: ( color online) Structure of the Si02 nanorod (left), and quantum stress - strain curve 
showing bulk solid properties (right) 

is found to hold in this welL 

II. THE IDEALIZED QUANTUM SOLID 

In this section the problem of multiscale modeling is introduced by first stating clearly 
what is the system being modeled. While much of this introductory matter is familiar, it 
provides the context for assessing all approximations involved in the final model material. 
First, the quantum description in terms of ions and electrons is given and the limitations 
for practical implementation are noted. Next, reduced self-consistent descriptions are given 
for the ions and electrons, each coupled to the other through their average charge density. 
In this form the classical limit for the ions can be taken, allowing the use of MD simulation 
methods for their dynamics. The equation for the electrons is simplified in a different 
direction, exploiting the very different time scales for electron and nuclear motion (the 
Born-Oppenheimer limit). This final description constitutes the idealized quantum solid for 
which the subsequent multiscale modeling scheme is proposed. 

At the fundamental level a simple atomic or molecular solid can be described in terms 
of Ni "ions" with charge number Zi for the corresponding atoms of species i and a set of 
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A^e electrons, with overall charge neutrahty NiZi — Ng). To introduce the various levels 
of description it is useful to start with the density operator D for the system as a whole. 
Properties of interest A are given by the expectation value 

(A) = Tr,,iDA, (1) 

where the trace is taken over all electron and ion degrees of freedom. The density operator 
obeys the Liouville - von Neumann equation 

dtD + ^[H,D]^0. (2) 

The Hamiltonian operator H is comprised of the Hamiltonians and for the isolated 
systems of ions and electrons, respectively, and their interaction C/jg 

H^Hi + H, + Uie, (3) 

t/.= /drdr'^iMMl. (4) 
J \r — r I 

Here pi (r) and (r) are the ion and electron charge density operators. This is the most 
fundamental level for a quantum description of the system. For a pure state, Disa projection 
onto a state ^ which is determined from the corresponding Schroedinger equation. For small 
systems, this can be solved numerically by quantum Monte Carlo methods. At this level, the 
approximations involved (e.g., electron nodal location) can be considered mild and under 
control. However, quantum Monte Carlo methods are restricted to relatively small systems 
and the above fundamental description is not a practical method for application to larger 
systems, particularly if repeated calculation is required to follow the dynamics. 

In many cases the properties of interest are functions only of the ion degrees of freedom 
(e.g. structure), A — > A^, or they are purely electronic (e.g., optical), A — > A^. Then 
a description is possible in terms the reduced density operators for the ions and for the 
electrons, Di and Dg , resulting from appropriate partial traces over all the electrons or ions, 
respectively 

(A) = Tr.AA, {Ae) = Tr,D,A,, (5) 
A = Tr^D, A = TnD. (6) 
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Their equations follow directly from (j2I) 

dtD, + '-[H,, A] + ^ (t/.A - Af//) = 0, (7) 
dtD, + '-[H,, A] + ^ (a A - AA^ = 0. (8) 

where the potential energy operators coupling the ion and electronic degrees of freedom are 
A = / dvdv'y^^p, (r) pe (r') , Pe (r) = (TrePe (r) D) {Tr,D)-' , (9) 

A = / drdr'y^^p, (r) pe (r') , p, (r) = (Tr.p^ (r) D) (TnD)-' (10) 

This is similar to the microscopic ion - electron coupling of Q except that now the electron 
charge density operator pe (r) is replaced by its conditional average, pe (r), in the equation 
for Di, and the ion charge density pi (r) is replaced by its conditional average, pi (r) in the 
equation for D^. The description ((7j) and (jS)) is still exact but formal since these average 
charge densities are not determined by these equations alone. The simplest realistic approx- 
imation (mean field) is to neglect the direct correlations in the charge densities pi (r) and 
Pe (r), i.e. replace D — > A A in ((7j) and (fTUj) to get 

Pi (r) TriPi (r) A = Pi (r) , pe (r) ^ Tr^pe (r) A = Pe (r) • (H) 

As a consequence, the potentials A and A become Hermitian and the average charge den- 
sities are now self-consistently determined by ((Zj) and (jHl). Self-consistency is required since 
De and A are functionals of the average charge densities p^ (r) and p^ (r), respectively. 

The advantage of this reduced description is that the ions and electrons are described by 
separate equations, reducing the degree of difficulty of each. More importantly, this separa- 
tion allows the introduction of appropriate approximations for each. The large differences 
in electron and ion masses imply corresponding differences in time scales and thermal de 
Broglie wavelengths. Consequently, the equation for the ions admits a classical limit for the 
conditions of interest, while that for the electrons does not. The classical limit of equation 
((Zj) becomes 

9iA + {(A + A[Pe]),A} = 0, (12) 

where {, } now denotes a Poisson bracket operation, and A is a function of the ion positions 
and momenta, A Di ({RjQ,}, {Pia}, t) (here a denotes a specific ion). This equation now 
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can be solved accurately and efficiently by molecular dynamics simulation methods, even for 
large systems. This is an essential step in almost all practical descriptions of bulk materials 
whose importance cannot be overstated. 

Implementation of still requires calculation of the potential energy Ui [pg], which in 
turn requires determination of the electronic charge density from De- However, the general 
solution to (jHI) is a formidable problem: determination of the dynamics of Nf, electrons self- 
consistently in the presence of a changing ion charge density. Two simplifications are made 
to bring this problem under control. First, it is recognized that r (dtDi) D^^ << 1, where r 
is the time scale for changes in the electron distribution function De. Second, it is assumed 
that only the lowest energy state contributes to at any given time. The first approxi- 
mation constitutes the Born-Oppenheimer approximation while the second approximation 
is appropriate for most structural studies, but must be relaxed for optical studies involving 
electronic transitions. In principle, (fT^ is solved in time steps At. At each time step the 
electronic charge density is computed for the ion configurations at that time step in order 
to recompute new forces for the next time step. The electronic charge density calculation is 
a ground state eigenvalue problem for the given ion configuration. 

To summarize, the final description of this idealized solid, it consists of a set of 
point ions governed by the classical equation ()12|) for their probability distribution Di — >■ 

A({Ria},{P^a},t) 

9iA + {(i/. + f/aPe]),A} = 0, (13) 

Ui = j (irdr'|— (r) (r', t) , p^ (r, t) = Tr^pe (r') A [Pi] (14) 
The average electron density p^ (r) is determined from the ground state solution to (jHI) 

+ A,A] =0, (15) 

Ue = J drdr'j-^^p, (r, t) p, (r') , p, (r, t) = Tr,p, (r) A [Pe] (16) 

The analysis proceeds stepwise. The classical equations (jl3|) are solved analytically in dis- 
crete time steps for the atomic coordinates and momenta. At each step the electron problem 
()15|) is solved for the electron ground state using the ion configuration at the previous time 
step. From this ground state the electron charge density pe (r) is determined. This gives the 
potential energy function Ui and consequently the forces required to change the ion positions 
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and momenta at the next time step. The process is repeated with a new electron charge den- 
sity calculated at each time step using the new ionic configurations. All electron correlations 
are accounted for quantum mechanically in the eigenvalue problem; all atomic correlations 
are determined classically through direct solution of Newton's equations. All structural 
properties of interest can be calculated since the phase points for the ions are known at all 
times. The dynamics of the electrons is only coarse-grained as they are " slaved" to the time 
dependence of the ions. 

The above defines a "quantum molecular dynamics" representation of the idealized 
quantum solid. The semi-classical approximation for the ions, and ground state Born- 
Oppenhiemer approximation for the electrons are relatively weak under most conditions of 
interest. The resulting description allows an accurate classical treatment of the ionic struc- 
ture while retaining relevant quantum chemistry for the interatomic forces due to electrons. 
If it could be implemented in practice for bulk systems of interest over reasonable time in- 
tervals there would be little need for multiscale modeling. With MD simulation the solution 
to (fTT^jl once Ui has been provided is straightforward. So the problem has been reduced to 
a determination of the electron charge density. Unfortunately, the solution to using 
realistic quantum chemical methods for even a few hundred ions at each time step becomes 
prohibitively time intensive. 

III. THE FORMAL PARTITION AND COMPOSITE SOLID 
A. The Representative Classical Solid 

In many cases of interest (e.g., equilibrium structure, thermodynamics) the computa- 
tional limitations of quantum chemical methods can be avoided through a purely classical 
representation of the solid, avoiding the intensive electron charge density calculation at each 
time step. This entails an idealization that has many variants. It consists of the representa- 
tion of the true potential energy function Ui ({RjQ,}) in (fT^ by a suitably chosen function 
Uc ({RjQ,}). Consequently, its form does not need to be computed at each time step and the 
speed and efficiency of classical molecular dynamics is not compromised. 

The problem with this approach lies in the choice for Uc ({RjQ,}). In principle, an exact 
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mapping for some fundamental property such as the free energy F can be imposed 

F[Ui] ^ F,[U,]. 



(17) 



where Fc is the corresponding classical functional. Generally, such exact methods can be 
inverted only perturbatively and lead to a sequence of effective many ion interactions involv- 
ing increasingly more particles. A more practical method is to assume pairwise additivity 
for effective point "atoms" 

Uc ({Ria}) ^ ^ E E E ^fci (\^ka - R,-/3|) , (18) 

where k, i label the species (ion or electron). The exact determination of the pair potentials 
Vkj (|Rjta ~ R-j^l) is now much more restrictive as not all properties of interest have such a 
representation. Nevertheless pair properties such as the radial distribution functions gkj{\r\) 
might be used to determine the pair potentials. As the gkj {\r\) are unknown and difficult 
to calculate, the inversion is again difficult and not practical. Instead, experimental data 
is often used to fit a parameterized functional form chosen for the pair potentials. At this 
stage control over the approximation is lost and the method becomes phenomenological. 

This phenomenological approach has been and remains a valuable tool of materials sci- 
ences. However, in the context of multiscale modeling it must be reconsidered carefully. 
First, it is recognized that there are local domains far from equilibrium where a purely clas- 
sical potential cannot apply because of the inherent quantum chemistry active there (e.g., 
charge transfer and exchange). Conversely, there are large complementary domains in near 
equilibrium states where representation by an appropriate classical potential is possible. 
Multiscale modeling constructs distinct classical and quantum models of these subsystems 
and then requires fidelity at their interface. This is a severe test of the modeling assump- 
tions in each subsystem. The approach proposed here confronts this issue directly in the 
construction of appropriate pair potentials for the problem considered. 

B. The Composite Quantum/Classical Solid 

It is presumed that there is some method for identifying domains within the solid for which 
quantum chemical effects should be treated in detail. The quantum solid is then represented 
as a composite of two domains, the larger bulk in which a classical representation is used 
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and a smaller "reactive" domain in which the original quantum description is retained. The 
objective of the modeling described here is therefore to construct the potential function Ui 
such that it gives an accurate description in both the reactive and non-reactive domains. 
This has two components, the determination of a pair potential for the forces on ions in the 
classical domain, and an accurate calculation of the charge density for the forces on ions in 
the quantum domain. 

The classical and quantum domains are defined by labelling all ions as either classical 
or quantum, and associating spatial domains with the coordinates of each, denoted {Rcio} 
and {Hqia} respectively. The quantum domains are assumed to be small, to allow practical 
calculation of the electronic structure. In principle there could be several disconnected 
quantum domains, but to simplify the discussion we consider only one. It is assumed that 
initially the two sets are contiguous with a smooth interface and that diffusion or migration 
between them is not significant over the times of interest. In addition to the ions in the 
quantum domain, there are m electrons where m is determined by a condition on the charge 
of the quantum domain, here taken to be neutral. The total average electron charge density 
is then decomposed as 

p, (r, t) = (r, t) {XQ (r) + xc (r)) = (r, t) + p,, (r, t) , (19) 

where xq (i") and xc (r) are characteristic functions for Q and C. The boundaries of the 
quantum spatial domain Q are constrained by the choice of ions and total electron charge 

"T-e = ^ drp^g (r, t) , (20) 

where Q encloses all {Rgja} (this leaves some flexibility to define smooth interfaces). The 
complement of this is the classical domain C. This gives a corresponding decomposition of 
the total ion potential energy function for the ions Vi (ion - ion Coulomb interactions plus 
Ui of (HH)) so that (O becomes 

dtD, + {{K, + V,, + V,,) , A} = 0, (21) 

where Ki is the ion kinetic energy and 

V.c = U dvp, (r) f / rfr'^-L- (p, (rO - 5 (r - r') + Pec (r', t)) + j dr' (p, (r') + p,^ (r', t)) 

(22) 
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and 

V., = U dvp, (r) [ f rfr'^-L- (p, (r') - 6 {r - r') + p,^ (r, t)) + f dr'j-^ (p, (r') + p,, (r', t)) 

(23) 

The first terms of tlie integrands on tlie riglit sides of and fl2H|l represent the interactions 
of the ions with a given subsystem in the presence of the average electronic charge density 
of that subsystem. The second terms represent the interactions of those ions with their 
complementary subsystem. Half of the ion - ion potential energy between the two subsystems 
has been associated with each potential in this decomposition so that the total force acting 
on the quantum domain by the classical domain is equal and opposite to that on the classical 
domain due to the quantum domain. 

The potential Vic is due to ions. By definition these ions are in near equilibrium states 
and therefore this part of the potential should be represented well by classical pair potentials 
of the form (jlSj) . with appropriately chosen parameters as discussed below. The potential 
Viq is due to ions in the quantum domain. As this is the domain that can be far from 
equilibrium the average electronic charge densities must be calculated in detail from the 
quantum description (fTK|) . There are two parts to this charge density affecting the ions in 
the quantum domain, that due to the m electrons of the quantum domain p^^{r,t), and 
that due to the surrounding classical domain pj (r) + p^^ (r, t). The equation governing the 
m electrons of the quantum domain is coupled to this same classical domain average charge 
density. The scheme proposed here consists of modeling this charge density Pg^ (r, t) as an 
accurate and practical representation for the environment of the quantum domain, allowing 
calculation of p^^ (r, t) and therefore determining both Viq and Vic- 

In summary, the multiscale model is obtained by replacing the ion - ion contribution in 
the classical domain (the first term on the right side of by suitably parameterized pair 
potentials, and constructing an accurate approximate calculation of the electron density in 
the quantum domain. The latter entails a representation of the effects of the electron density 
in the classical domain on charges in the quantum domain. In this way the potentials of 
()22|) and (jSBI) are entirely determined. The details of this construction are addressed in the 
following sections. 
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IV. DESCRIPTION OF THE CLASSICAL SUBSYSTEM 



A. Construction of the pair potentials 

The proposed method for constructing a classical pair potential for use in multiscale 
modehng has several components: 1) it should be constructed for accuracy of the specific 
properties to be studied, 2) it should be "trained" on quantum data generated in the same 
way as for the quantum domain at its interface, 3) it should predict accurate equilibrium 
structure, but include training on appropriate near equilibrium states as well, and 4) sim- 
plicity of form for parameterization and implementation in MD codes should be maintained. 

The steps in constructing a potential are the following. First the specific quantum method 
to be used in the multiscale modeling is identified (below, transfer Hamiltonian or density 
functional) and applied to a large cluster or representative sample of the solid to be modeled. 
The forces on atoms for both equilibrium and near equilibrium states are then calculated 
quantum mechanically. Next, a simple functional form for the pair potentials with the 
correct physical shape (e.g., one of the existing phenomcnological forms) is chosen and 
the parameters controlling that shape selected for optimization. The forces on the ions 
are calculated for these chosen potentials at the configurations used in the quantum force 
calculations, and compared with those quantum forces. The parameters are adjusted for 
a best fit to the quantum force data. When a good fit has been obtained, the potential 
energy at equilibrium is tested for stability using a gradient algorithm to establish that a 
local minimum of the potential has been obtained. Finally, the property of interest (e.g., 
some linear response) is tested by comparing its calculation using MD simulation with the 
fitted potentials and that with the quantum forces. If necessary, the fitting procedure 
can be repeated with differing weights for the quantum force data in equilibrium and near 
equilibrium states, or other additional input from the quantum calculations. 

A more detailed discussion of the flexibility in this procedure is described for the example 
of the next section. It entails some art as well as science in choosing the optimization methods 
and applying them. In that example, it was found that some repeated combination of genetic 
algorithms and scaling provided the accuracy required. The primary result of this approach 
is a classical pair potential which gives the properties of interest by design, that match those 
being calculated quantum mechanically across the interface - for near equilibrium states. 
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B. Effects of the environment 

The above construction describes the modehng of the first term on the right side of ()22|) 
where both the ions and the average electron charge density refer to the classical domain, 
assumed in a near equilibrium state. The second term involves a coupling to ions and average 
charge density (now the ground state charge density) in the quantum domain (r, t). This 
charge density is computed for the forces in the quantum domain (next section) and hence the 
second term of (j^^ is known as well. Thus, the potential energy for the ions of the classical 
domain is determined from synthesized pair potentials among the ions and electrons of the 
classical domain, plus an interaction with the electrostatic potential of the ions and average 
electron charge density of the quantum domain. However, in the examples discussed below 
the coupling of the classical domain ions to the quantum domain is simplified by using the 
same pair potentials as for ions within the classical domain. This is expected to be quite 
accurate if the quantum domain charge density is not distorted very much from the near 
equilibrium states, since the potentials are trained to be equivalent to the charge density in 
the near equilibrium domain. 

V. DESCRIPTION OF THE QUANTUM SUBSYSTEM 

The quantum domain is a subsystem of m electrons localized about the designated ions 
defining that domain. Consider the reduced density operator for m electrons defined by 

Jjim) ^ Tr^e-mj^^ ^24) 

More specifically, this partial trace is defined in coordinate representation by 

(ri, ..,rm| D^J^^ \r[, ..,r^) = J drm+i.-dvN, (ri, ..,rm,r„, ..VN,\De \r\, ..,r'^,rm, ..tn,) (25) 

Clearly the full exchange symmetry among all A'^e electrons is preserved. However, this 
reduced density operator is not specific to the quantum domain defined above only. For 
example, the diagonal elements (ri, .., r,m| -Dg |ri, ..,rm) give the probability density to find 
m electrons at the specified positions, and the latter can be chosen anywhere inside the 
system. Thus, only when the positions are restricted to the quantum domain Q does this 
reduced density operator represent the electrons of that domain. Similarly, if p^*") (r) is the 
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charge density operator for m electrons its average is 

^ r^^jrn) jj^ ^ Tr^Tp'r^ (r) D^J^^ (26) 

where the trace in the second equahty is over m degrees of freedom. This average 
charge density represents the average contribution of m electrons at any point r. Both 
{ri, ..,Vm\ Dq\ri, .., Ym) and p^'") (r) change with r since there is an absolute reference back- 
ground set by the functional dependence on the ion charge density. Consequently, in all of 
the following discussion of this section (ri, .., rm| |ri, ..,Vm) and p^™^ (r) are considered 
only for positions within the chosen quantum domain. Accordingly, this coordinate rep- 
resentation {|ri,..,rm)} defines an m electron Hilbert space of functions defined over the 
quantum domain Q. In this context D^^^ becomes the reduced density operator for the 
quantum domain and p^™'-' (r) its charge density operator 

D(-) ^ D„ p(-) (r) pe, (r) . (27) 

The equation determining the reduced density operator for the quantum subsystem fol- 
lows directly from this definition and Eq. (fTH|l for 

'-[{K, + Ve,),D,] = Q. (28) 

where Ke is the kinetic energy for the m electrons and 

Ve, = \l drp,, (r) ( I rfr'^— (p,, (r') - J (r - r') + p, (r')) + / dv' (p,, (r') + p, (r')) 
Z JQ \Jq |r — r I Jc |r — r I 

(29) 

(The contribution from p^^ (r') in the second term actually should be a conditional average; 
the same mean field approximation as in (jlip has been introduced here for consistency). The 
first term on the left side of (j29|l describes the isolated quantum subsystem, composed of 
the Coulomb interactions among the m electrons and their coupling to the average charge 
density of the ions in the quantum domain. The second term is the interaction of these 
electrons with their environment, the total average charge density of the classical domain. 

There are two distinct types of contributions from this charge density of the classical 
domain. The first is associated with a subset of ions at the border of the QM/CM domains 
which describe chemical bonds in the full quantum solid. These ions locate regions where 
there is a highly localized electron charge density shared with the quantum domain, including 
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both strong correlation and exchange effects, locahzed and non-uniform. The second type 
of contribution is associated with the remaining border ions and those charges more distant. 
In this second case the quantum correlation and exchange effects of the first type are much 
weaker, and the dominant effect is that of a polarized charge density due to the ions and 
electrons. 

These two types of effects of the environment can be identified by separating the total 
charge density for the classical domain in (j^^ into domains centered on the border ions 
responsible for bonding, and their complement 

(r) ^ p,, (r) + (r) = 5: X (|r - Ral) Pe (r) + Ap, (r) (30) 

where it is understood that r is in the classical domain. The set of border ions for which a 
bond has been broken in identifying the quantum domain is denoted by V. Also, x{\^ ~ R-al) 
is a characteristic function specifying a domain centered on Rq, such that it does not overlap 
neighboring ions. Its size is taken large enough to incorporate the bound electrons forming 
the " atom" for this ion in the quantum solid. The second term Ap^ (r) is the charge density 
for all remaining ions and electrons of the classical environment. 

A. Coulomb effects of environment 

By its definition, the averaged charge density Ap^ (r) does not include the contributions to 
chemical bonding with the quantum domain required for its valency. Hence the electrostatic 
potential associated with it can be expected to have a regular multipole expansion 

/ dr'- -Ap, + .. 31 

JQ |r — r I 

The leading monopole term is zero due to charge neutrality of the classical subsystem, and 
the dipole moment d for the entire environment is 

d=J dr'r'Ap^ (r) . (32) 

These results are still quite formal but they provide the basis for the phenomenology 
proposed for this part of the modeling: Replace all effects of the classical environment on the 
quantum system, exclusive of bonding, by an effective dipole representing the polarization 
of the medium by the quantum domain. The origin of this dipole in the above analysis 
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shows that in general it will depend on the geometry and the state of both the classical and 
quantum domain (e.g., it will change under conditions of strain). In the phenomenological 
application of this prescription the dipole must be supplied by some simpler means since the 
electronic contribution to Ap^ (r) is not known (for example, see referenc^). 

B. Electron exchange with environment 

The contributions from (r) in the regions where bonds have been cut require a more 
detailed treatment. Clearly, a necessary condition is that valence saturation should be 
restored in the quantum domain. However, that is not sufficient to assure the correct charge 
density there nor the correct forces within that domain. Instead, the charge density near 
the border ions responsible for bonding should induce a realistic charge density within the 
quantum domain. To see how this can be done first write the contribution from one such 
border ion to (jHUj) as 



The potential (ppa (r) represents the actual electrostatic and exchange effects of the nucleus 
of the border of the quantum domain. This is comprised largely of the ion at the site plus 
its closed shell electrons distorted by exchange and correlation effects of the site in the 
quantum domain with which it is bonding. Consequently, an appropriate pseudo-potential 
is introduced at each such border ion whose behavior in the direction of the quantum system 
is the same as that of the actual ion. These ions plus their pseudo-potentials will be called 
"pseudo atoms". 

To accomplish this, an appropriate candidate for the pseudo atom is chosen based on 
valency of the particular pair of border ion and its neighbor in the quantum domain. Next, 
a large molecule or small cluster containing that bonding pair is chosen for training the 
pseudo atom. The bond is then broken and the relevant member of the pair replaced by the 
pseudo atom. The training consists of parameterizing the pseudo atoms' effective potential 
to give the same forces as in the original cluster. For example, the pseudo atom might 
consist of an ion plus its closed shell electrons chosen to satisfy the valency for the bond 
broken. The adjustable parameters could refer to a characterization of the closed shell 
electron distribution. In this way, it is assured that the pseudo atom not only gives the 




(33) 
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correct saturation of the dangling bond but also reproduces the forces within the cluster 
and hence gives a realistic representation of the charge distribution between the chosen pair. 
For the small training cluster (e.g., Figure IH)), forces are determined both at equilibrium 
and under strain to within about one tenth percent. The primary assumption is that the 
bonding of interest is a local effect, so that the pseudo atom trained in the cluster will have 
a similar accuracy when used in the bulk solid . The training depends on the particular 
quantum method used to describe the solid, and is illustrated in more detail for the specific 
example considered in the later sections. 

In summary, the environmental effects on the quantum system are accounted for approx- 
imately by a dipole representing its polarization and a collection of pseudo atoms located 
at the sites of ions where bonds have been cut 

Veg^U drp^, (r) ( [ rfr'^- (p., (r') - 5 (r - r') + p. (r')) + ^ + E (r)) • 

(34) 

This model now allows solution to ()28p for the electron distribution the quantum domain, 
including its coupling to the classical environment. Then Peg (r) is calculated from (j^UI) and 
(f?fjl . Finally, the desired potentials Vic and Viq of and are fully determined 

1 I r /■ 1 / N 

= ^ E E E - + / ^^^^ / ^^'^^ ^^') + "P-'i ^)) ' (35) 

^ i,j a [3 •'^ JQ \r r \ 

and 

V^g = l [ drp, (r) ( / rfr (p,; (r') - 5 (r - r') + p,, (r, t)) + ^+J2 <Ppa (r) ) • 

^ JQ \JQ |r r I r ) 

(36) 

The solution to the classical ion motion ()2H1 can then proceed via MD simulation. This 
completes the proposed scheme for multiscale modeling of the idealized quantum solid. 
Rather than critique further the assumptions made in this abstract context, the scheme is 
discussed in more detail for a specific application to an 5*^02 nanorod. 

VI. 5^02 NANOROD - A CRITICAL TEST 

The modeling scheme of the previous sections is now illustrated in detail and tested 
critically for a model solid consisting of a silica nanorod^^^ shown in Figure ^ It consists 
of 108 Si and O ions with the stoichiometric ratio of one silicon to two oxygens, a stack 

18 



of six SiqOq rings sharing both above and below a ring of oxygen atoms. To saturate 
overall valency the nanorod is terminated with end cap rings whereby each silicon atom 
of the terminating ring is connected by bridging oxygen or two interstitial oxygens. The 
size of this nanorod can be readily adjusted by adding or removing Si^O^ planar rings with 
corresponding O atoms. Many of the results reported here have been studied for larger rods 
as welP^. As noted in the Introduction, this system exhibits strain to fracture (Figure 1), yet 
it is small enough for practical application of the chosen quantum mechanics to the whole 
rod. Consequently, the results of modeling it as a composite rod can be tested quantitatively 
against the "exact" results. That is the objective of this section. 

The decomposition of the rod into classical and quantum domains is accomplished by 
identifying the ions of one of the rings near the center as the QM domain. The rest of the 
ions define the classical domain. The actual geometry of the interface is chosen by charge 
neutrality of the quantum domain, as given by Eq.fpUj). 

A. Transfer Hamiltonian Quantum Mechanics 

At the quantum level, each possible choice for the quantum chemical method entails 
different approximations used to optimize the accuracy and speed of the calculations. Each 
approximation has its own advantages and limitations, but a general embedding scheme 
should be insensitive to these. Much of the previous work on multiscale modeling is based 
on tight binding models^^ because of their ability to treat hundreds of atoms in the QM 
domain. However, the Hamiltonian in these models is oversimplified and cannot account for 
charge transfer, which is essential to bond breaking. A more complete quantum description 
is considered here. It is appropriate to stress at this point that the objective of the following 
sections is to demonstrate and test the modeling scheme for a given quantum description. 
That is a separate issue from the question of how accurate that description may be. To test 
how robust the proposed embedding scheme is, two different methods are chosen for the 
underlying quantum mechanical description. The first method is the Transfer Hamiltonian 
(TH) - Neglect of Diatomic Differential Overlap (NDDO)^. The second method is the 
Born-Oppenheimer local spin density^^ within density functional theory using a generalized 
gradient approximation. The results obtained from the TH are described first. 

The goal of the TH strategy is to provide forces for realistic MD simulations of the quality 
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of the coupled cluster singles and doubles method^S in computationally accessible times^S. 
The TH is a low rank, self - consistent, quantum mechanical single particle operator whose 
matrix elements are given in terms of parameterizable functions. The functional form of the 
TH can be chosen to be of any semi-empirical functional form, but here the NDDO form is 
chosen because of its better qualitative description over Intermediate Neglect of Differential 
Overlap (INDO) Hamiltonians^. The NDDO approximation restricts charge distribution to 
basis functions on the same center, so this has at most two-atom interactions. Because of its 
NDDO form it is anticipated that the parameters fit to coupled cluster data are saturated 
for small cluster size and that these parameters can then be transferred to extended systems. 
This transfer might be considered an example of serial multiscale-modeling. 

To apply the TH to silica systems for studying complex processes such as fracture or 
hydrolysis where the role of electrons is important, Taylor et al.'^'^ parameterized the TH for 
the following reactions Si2HQ — > 2SiH^, Si20'rHQ SiiOH)^ + SiOiOH)^ (radical mode 
of dissociation), and Si207HQ — > ^ Si{OH)^ +^ OSi{OH)^ (ionic mode of dissociation). By 
describing the two different dissociation pathways for Si20iHQ, electron state specificity is 
systematically introduced into the TH model. Since fracture is expected to involve a variety 
of dissociations this TH model is expected to be appropriate for study of strained conditions 
as is done here. Thus the parameterized TH model is expected to provide accurate results for 
such processes, while requiring only the computational intensity of semi-empirical methods 
which are orders of magnitude less demanding than coupled cluster calculations on the same 
systems. In the next few subsections, the charge densities and forces that arise from the TH 
are the reference data by which the proposed method for multiscale modeling is judged. 

B. Pair Potentials for the Classical Domain 

The method that we propose for constructing a classical pair potential has been previously 
presented in detail^, so only a brief summary is given here. The functional form for the 
potential is chosen to be one that is commonly used in the literature. In the case of silica, 
this is the TTAM^^ and BKS^ forms which are pair potentials comprised of a Coulombic 
interaction, a short-range exponential repulsion, and a van der Waals attraction (the last 
two terms are collectively known as the 'Buckingham' potential) 




(37) 
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There are 10 free parameters as the charge neutrahty of the Si02 unit requires that qsi = 
—2qo, so that aij, bij, and Cij for each atom pair and one charge define the potential. For 
the purposes here, this potential is used to generate the dynamics (MD simulations) so 
that forces have been used as the property "training" the potential (determining the free 
parameters). More specifically, force data from the silica nanorod calculated with the TH 
at equilibrium and a few small strain configurations have been used. 

Fitting a set of ten parameters from hundreds of data points requires an optimization 
approach that is capable of exploring a large parameter space efficiently without being 
trapped in local minima. For this task, a genetic algorithm (GA) has been chosen as detailed 
in reference^^. The GA generates a number of sets of parameters of a given fitness which are 
then tested using a standard geometry optimization technique (BFGS'^^) to ensure that the 
parameterization gives a stable equilibrium configuration. The parameter set that gives a 
geometry closest to the TH nanorod is then chosen and rescaled to reproduce the equilibrium 
structure to within a few percent. (See the appendix of reference^^ for details of the rescaling 
used.) Table IJ presents the final parameters obtained compared to the TTAM and BKS 
parameters. Table ITTl shows that the new classical potential reproduces the nanorod structure 
to within a few percent compared to larger discrepancies with the two standard potentials. 
The shapes of these three potentials are shown in Figure |21 for the Si-0 pair interaction. 

Finally, the small strain elastic behavior of this new potential is studied. This is accom- 
plished using MD simulation of uniaxial strain along the long axis (z) of the silica nanorod 
at a strain rate of 25 m/s. The integration is done using 2 fs time steps and a temperature 
of 10 K enforced by velocity rescaling. Figure El shows the resulting stress-strain curves for 
small strains up to 5%. It is noted that the new potential indeed reproduces the linear 
elastic response (Young's modulus) behavior of the underlying TH while the TTAM and 
BKS potentials are somewhat stiffer. 

A key to successful multiscale modeling lies in the consistent embedding of a QM domain 
in its classical MD region. This requires that the CM region have the same structure and 
elastic properties as the QM domain. The potential obtained in this section meets these 
criteria for the silica nanorod. Essentially, a classical representation of the TH silica nanorod 
has been found for both equilibrium and near-equilibrium configurations. 
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TABLE I: Potential Parameters. 
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FIG. 2: (color online) Si-0 pairwise interactions for the three potentials 
C. The Quantum Domain 

A proper description of the quantum domain requires incorporating two kinds of envi- 
ronmental effects: a short-range electronic exchange interaction and long-range Coulombic 
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TABLE II: Structure of the Nanorod from the Different Potentials. 
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FIG. 3: (color online) Stress-strain curves for the different potentials and for TH quantum mechanics 

interactions. Most previous studies of QM/CM simulations have focused on the former and 
have neglected the long-range interactions. It will be shown here that both kinds of interac- 
tions must be taken into account for an accurate description of forces and charge densities 
in the QM domain. This subsection first reviews briefly the various termination schemes 
used in the literature to treat the bond cutting region. Then, the method proposed here 
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for construction of a pseudo-atom to saturate dangling bonds is described, together with 
modehng the rest of the CM environment by lowest order multipoles. The charge densities 
and forces in the QM domain obtained from our scheme of pseudo-atoms and dipoles with 
those obtained from the conventionally used link atoms (bond saturation with hydrogen 
atoms) and dipoles are compared to the "exact" TH reference data for the entire nanorod. 

Over the past decade there have been numerous proposals for different types of termi- 
nation schemes to accommodate dangling bonds. Only a few of the more relevant ones are 
noted for context, (i) Link Atom method: This commonly used method is the Link Atom 
method, presented by Singh and KoUman^S. In this method, hydrogen atoms are added to 
the CM side of broken covalent bond to satisfy the valency of the QM system. There are 
many variations within the implementation of the LA method, for example, the double Link 
Atom method^°, the Add- Remove Link Atom method^^ or the Scaled- Position-Link- Atom 
Method (SPLAM)^. In the present work, the hydrogen atoms are placed at a fixed distance 
of 0.97 A from the cut Si bond, (ii) Connection Atom method^: A connection atom is 
developed to saturate a C-C bond, such that the connection atom mimics the effect of a 
methyl group. This connection atom interacts with the QM atom quantum mechanically 
and the interactions with the CM atoms are handled classically using a carbon force field. 
The parameters of the connection atom are determined using semi-empirical methodsAi such 
as AMI, MNDO or PM3 designed to reproduce theoretical QM data for energies, geometry 
and net charges. About 30 different methyl hydrocarbons were used as reference molecules. 
The parameters adjusted are the orbital exponent , one-center one-electron energy , one 
center two-electron integral , resonance parameter , and repulsion term . The mathematical 
functional forms of these parameters can be found in any book on semi-empirical theory^ 
or papers by Dewar and Theil^. 

There are numerous other termination schemes for the QM/CM boundary like the 'pseu- 
dobond' scheme^, TMOMM'^^i,^ and 'ONIOM'iSa,^^!,^ procedures , 'effective group potential 
(EGP)'-^. However these methods are not discussed here. 

The method proposed here will be referred to as the "pseudo-atom" method. In the 
present case it is based on the TH approach for saturating a bond terminating in Si (in 
silica systems). The training cluster is a pyrosilicic acid molecule (Figure 0]). The part 
of the molecule within the dotted lines is replaced by a fluorine (F) atom whose NDDO 
parameters are then adjusted to give the correct QM forces (which implies correct geometry 
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as well) and charge density in other parts of the molecule (outside the dotted lines). The 
pseudo atom is placed at the same position as the neighboring CM atom (O in this case) in 
the bond being cut and hence geometry optimization as in the LA method is not required. 
The specific NDDO parameters modified are: one-center-one-electron integrals , Coulomb 
integrals, exchange integral, and two-center-one-electron resonance integrals^. This method 
is similar to the previously described connection atom developed in reference^. However 
the advantage of using pseudo-atoms is that they are trained to give the correct QM forces 
and charge densities for both the equilibrium Si-0 bond being cut and for small distortions 
(up to 3 — 4% from the equilibrium) in the Si-0 bond length as well. This allows use of 
the pseudo-atoms even while studying dynamics in the system. To test this method, the 




FIG. 4: (color online) Training of Pseudo-atom on Pyrosilicic acid 

actual charge density from the TH method is calculated instead of the more commonly used 
MuUiken populations^, which are known to have several common problems ( e.g., equal 
apportioning of electrons between pairs of atoms, even if their electronegetavities are very 
different). 

The remainder of the environment is represented as two dipoles for the top and bottom 
portions of the rod (Figure EJ. The values of the dipole have been calculated using the 
TH-NDDO charge density for these two portions of the rod. These domains are taken to be 
charge neutral, but are polarized by the presence of the QM domain. The validity of the 
approximation by dipoles has been checked by comparing the force on an Si nuclei of the 
ring due to all charges and that due to the dipole, with excellent agreement. Constructing 
the full system as a composite of a QM domain terminated by pseudo atoms and embedded 
in dipolar fields leads to the forces displayed in Figure IHl It is found that the normalized 
difference of charge density is reproduced to within 0.1% in the plane of the ring^i. 

This scheme was found to be applicable to both equilibrium and strained configurations, 
as illustrated in Figure IHl for the following CcLSGSi cL ) equilibrium, b) the ring of the QM 
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FIG. 5: (color online) Embedding Scheme: Approximating the CM region by pseudo-atoms and 
dipoles 

domain radially expanded by 5%, and c) a distorted ring in which one Si atom is radially 
pushed out and one Si pushed in. Also shown are the corresponding results for a longer 10 
ring rod (cases d),e), and f)). 




FIG. 6: Forces on Si nuclei for various cases studied 

The pseudo-atoms and the dipoles reproduce the force in the QM domain within 1% 
error. In contrast, the LA placed at 0.97 A from the terminating Si atom and along the 
Si-0 bond lead to a large force on the Si atom at equilibrium, and generally poor results 
for the strained cases. One might argue that the LA method could be improved if placed at 
"optimal" positions. To test this, both the Si-H bond length as well as the alignment of the 
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LA was varied to give minimum force on the Si atom. This was obtained when the LA is 
placed at a distance of 1.45 A from the sihcon and the bond angle is decreased by about 5 
degrees. Although the LA at this position gives forces comparable to those for pseudo-atoms 
plus the dipoles, it fails to reproduce the correct charge densities. 

D. The Composite Rod 

The composite rod is built by embedding the QM region in its CM environment as 
described in Section 111. The forces on the atoms in the CM region are calculated from the 
pair potential developed in IV B above, while the forces on the atoms in the QM domain 
are obtained from the charge density calculated by the TH method as described in VI C . 
An important test of the composite rod is its indistinguishability from the TH rod for near 
equilibrium states, i.e. its structure and elastic properties. 

The structural properties (bond lengths, bond angles) have an accuracy comparable to 
that of the purely classical rod described using this pair potential (Table ITT|) . Consequently, 
attention here will be focused on the elastic properties. Figure [7| shows the stress-strain 
behavior of the nanorod obtained from three different methods: (i) quantum mechanics TH 
method for the entire rod, (ii) pair potentials for the entire rod, and (iii) the composite rod 
constructed as described above. The three overlaying curves (measured at O.OIK) indicate 
that the composite rod is identical to the rod obtained from TH and the pair potential 
nanorod in terms of small strain elastic properties and structure. The stress-strain results 
shows the success of our multiscale method indicating that the composite rod is indistin- 
guishable from the underlying quantum mechanics for states near equilibrium. 

E. Notched Nanorod 

The elastic properties of the composite rod obtained from both of these potentials do 
not agree beyond 10% strain. This is because the pseudo-atoms, trained at regions only 
close to the equilibrium configuration, fail to give the correct charge densities at such high 
strains. This diagnosis was checked by comparing the charge density in the QM domain of 
the 12% strained rod with that of equilibrium configuration and a difference of 6% was found. 
Retraining of the pseudo-atoms improve the stress-strain performance of the composite rod 
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beyond a strain of 10% is possible, but this was not done because of the large strains involved. 
Real systems like glasses have many inherent defects which act as stress concentrators that 
cause the material to break at much lower strains than observed for the nanorod (with a 
yield stress of about 190 GPa). 
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FIG. 7: (color online) Stress-strain curve for the TH (quantum), new potential (classical), and 
composite rods 

To illustrate this, a defect notch was placed in the 108 atom nanorod by removal of an 
oxygen atom as shown in Figure |H1 The MD stress-strain curves for this notched rod were 
found using TH quantum mechanics, the trained classical potential, and the composite. We 
note that the presence of only a small defect can significantly reduce the yield stress of 
the material and make it more prone to fracture. The TH curve for the defect-free rod is 
plotted in the same figure to contrast the value of the yield stress. As can be seen there 
is a reduction of ~60 GPa in the yield stress. For the composite rod in this case, the QM 
domain was chosen to consist of 2 silica planes and the intermediate 5 oxygen atoms (see 
Figure IHI), so that the defect could be located in the QM region. The stress - strain curve for 
this composite notched rod agrees well with that for the TH quantum calculation up to 10% 
strain. Above 4% strain the curve for the composite notched nanorod follows that of the 
TH instead that of the trained classical potential nanorod, showing that the composite rod 
is representing the "real" material. This is exactly what is required of multiscale modeling. 
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FIG. 8: (color online) The notched nanorod and corresponding stress-strain curve 
F. DFT Quantum Mechanics 

A proper multiscale modeling procedure should be independent of the choice of underlying 
quantum mechanical method. To test the method proposed here, the analysis of subsections 
(B)-(D) is repeated using a density functional theory (DFT) instead of the Transfer Hamil- 
tonian as the quantum mechanical approximation. The primary results of this section are a 
confirmation that the isolation of the QM domain with pseudo-atoms and dipoles, plus the 
construction of a classical potential based on the DFT forces leads to accuracies of the same 
quality as those described already using the TH quantum mechanics. Hence, the multiscale 
modeling scheme is faithful to the chosen form for the underlying quantum mechanics in 
both cases. The DFT code used is a parallel multiscale program package, known as Born- 
Oppenheimer molecular dynamics "BOMD"^^. It is a generalized gradient approximation 
(GGA) within local spin DFT (LSDFT). A Troullier-Martin Pseudo-Potential^iS is used 
for the effect of the chemically inert core states on the valence states. The code employs the 
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dual space formalism for calculation of the DFT energy. A plane wave basis set and cut-off 
energy of 30.84 Rydbergs is chosen. 
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FIG. 9: (color online) Adiabatic stress curves for DFT, DFT-potential, and TTAM potential 

A pseudo-atom for the quantum domain is constructed based on parameterization of 
the Troullier-Martin (TM) pseudopotential, using a cut-off radius of 1.5 A. Once again the 
pyrosilicic acid molecule is chosen for parameterization of the fluorine-like pseudo atom and 
its position is constrained to be at the same place as the O atom (Figure |3)). Unlike the TH- 
NDDO method in which both the electron-ion and electron-electron interaction parameters 
could be changed , in DFT only the electron-ion interactions can be modified. The three 
options to alter these interactions for the F atom using the TM pseudopotential are: (i) 
the core charge on F, (ii) omission of the non-local part in the potential, and (iii) switching 
the local and non-local part between s and p orbitals. All three possible choices and their 
combinations were explored to find the optimal reproduction of forces on the terminated Si 
atom in the pyrosilicic acid. The best results were obtained when the core charge is 7.0 and 
the non-local part is omitted. This pseudo-atom was then applied to the composite nanorod 
constructed as above (Figure E)). The values of the dipoles were recalculated from the charge 
density obtained using the DFT results for the CM portions of the rod. The force on a Si 
atom in the QM domain was calculated for the rod in equilibrium and all the strained cases 
considered in subsection C. The results obtained from a DFT calculation on the whole rod 
are now taken to be the "exact" reference forces. It is found that the forces and the charge 
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densities in the QM domain can be generated using DFT pseudo-atoms and dipoles to 1% 
accuracy. 

Next, a new potential having the same form as TTAM is constructed to predict the same 
structure and elastic properties as for the QM rod. A GA with DFT force data up to 4% 
expansion followed by the scaling procedure is used as described in section VI B above to 
find the parameters for the potential. The charge on the ions is lower (as was found for 
TH quantum mechanics) than that given in the TTAM and BKS potentials. Also, the 
van der Waals interaction is much weaker than in those potentials. This is expected since 
the DFT forces fail to represent this effect, falling exponentially rather than algebraically 
with separation. The parameters for the DFT potential and a comparison of the resulting 
nanorod structure with that for the BKS and TTAM nanorods are given in reference^^. The 
agreement between the results from the DFT potential and DFT quantum mechanics is 
similar to that found in Table |n] for the TH quantum mechanics. 

A stress strain curve for the entire rod using MD with DFT quantum mechanics is com- 
putationally too intensive. Instead, only selected equilibrium and adiabatic strain configu- 
rations were calculated with DFT. The equilibrium structure was determined by sequential 
DFT calculations and nuclear relaxation to find the minimum energy configuration. The 
strained configurations were obtained from an affine transformation of the minimum energy 
configuration by 1, 2, 3, and 4 %, with a single DFT calculation of forces at each of the 
expanded configurations. Then the average force on the Si atom in each ring of the DFT 
nanorod was computed for these four cases. The stresses for adiabatic configurations using 
DFT for the entire rod, the constructed DFT potential, and the TTAM potential were then 
compared (the values of stress obtained using BKS potential are similar to those of TTAM 
potential). The results are shown in Figure IHl 

These results suggest that the multiscale modeling method proposed here is accurate for 
a wide range of possible choices for the underlying quantum mechanical method employed. 

VII. SUMMARY AND DISCUSSION 

The objectives here have been three-fold. The first was to describe the formal quantum 
structure for a real solid from which a practical model should be obtained through a sequence 
of well-identified (if not fully controlled) approximations. The second was to propose a 
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method for partitioning this structure into classical and quantum domains while preserving 
the properties of interest for the quantum structure in the replica composite solid. The final 
objective was a quantitative test of the proposed method by its application to a non trivial 
mesoscopic "solid", the Si02 nanorod. 

For the first objective, coupled Liouville - von Neumann equations for the reduced density 
operators for the ions and electrons of the system were considered. Each is coupled to the 
other through the mean charge density of the complementary subsystem. Determination 
of these charge densities then becomes the central problem for further analysis. For the 
heavy ion component analysis by classical MD simulation provides a practical and accurate 
approach. The necessity for a detailed quantum treatment of the electronic charge density 
is the primary bottleneck for progress. For bulk samples of interest, direct application 
accurate quantum chemical methods are precluded, so composite constructs with smaller 
quantum subdomains are a potentially fruitful compromise. Here, the isolation of a quantum 
subdomain was accomplished by defining the reduced density matrix for electrons in the 
vicinity of a selected small set of ions. The corresponding Liouville - von Neumann equation 
describes the dynamics of those electrons and ions coupled to the remainder of the solid 
(the classical domain). Its environment is entirely characterized by its mean charge density. 
This formulation and emphasis on the charge density is a guiding feature of the subsequent 
approximations and modeling. 

The specific method for constructing the composite sohd consists of modehng the charge 
densities of the environment for the two subsystems. For the environment of the quantum 
domain the charge density is separated into that border component responsible for electronic 
bonding across the border, and a longer range Coulomb interaction. The former is treated 
in detail by pseudo atoms trained on smaller clusters to provide highly accurate local forces 
and charge densities. The latter is represented by leading order terms in a multi-pole expan- 
sion. The interactions within the classical domain are described phenomenologically by pair 
potentials. These pair potentials are designed specifically to represent selected properties 
determined by the quantum chemical method assumed to represent the solid. This presumes 
conditions such that the chosen properties have such a classical representation (e.g., near 
equilibrium states), and the method constitutes a "tuning" of the potential to match the 
quantum description on the other side of the border. 

There are several critical tests of this approach, all requiring a benchmark system for 
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which the "exact" quantum data for the entire sohd is required. This system is provided 
here by the nanorod for which the global quantum calculations are possible. Then, choosing 
a central ring as the quantum subdoman and the remainder as the classical domain many 
quantiative tests are possible. Some of these provided here are: 

1. The charge density and forces in the central ring computed from the model are accurate 
to within one percent, both at equilibrium and strains up to 5%. 

2. The pair potential fit to quantum data was used to construct the entire nanorod. The 
resulting structure and elastic properties up to 5% strain were indistinguishable from 
those of the quantum nanorod. 

3. The proposed method was used to construct a composite classical / quantum nanorod. 
Again the structure and elastic properties were indistinguishable from those of the 
quantum nanorod. 

These and other results presented above constitute a demonstration that the multiscale 
modeling is not creating a new solid, but rather is faithful to the real system of interest. 
Similar (and perhaps more physical) results were obtained for the nanorod with a defect 
(missing Oxygen) . Finally, the entire multiscale analysis was repeated using a quite different 
choice for the underlying quantum mechanics with the same degree of accuracy. 

The predictions of multiscale modeling for conditions of interest, states far from equi- 
librium, rely on the assumption that the properties calculated in the quantum subdomain 
are indeed those of the given quantum structure. This means that the quantum imbedding 
and representation of the classical domain are "passive", i.e. no new physical effects have 
been introduced. The method described here manifestly satisfies this constraint for the 
benchmark nanorod, and provides some concrete support for its use in realistic applications. 
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